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ABSTRACT 



We report on the fourth phase of our study of slightly rotating accretion flows onto black 
holes. The main new element of this study is that we used fully three dimensional (3-D) numeri- 
cal simulations. We consider hydrodynamics of inviscid accretion flows. We assume a spherically 
symmetric density distribution at the outer boundary, but brake the flow symmetry by intro- 
ducing a small, latitude-dependent angular momentum. We also consider cases where angular 
momentum at large radii is latitude- and azimuth-dependent. For the latitude-dependent angu- 
lar momentum, 3-D simulations confirm axisymmetric results: the material that has too much 
angular momentum to be accreted forms a thick torus near the equator. Consequently, accretion 
proceeds only through the polar funnel, and the mass accretion rate through the funnel is con- 
strained by the size and shape of the torus, not by the outer conditions. In 3-D simulations, we 
found that the torus precesses, even for axisymmetric conditions at large radii. For the latitude 
and azimuth-dependent angular momentum, the non-rotating gas near the equator can also sig- 
nificantly affect the evolution of the rotating gas. In particular, it may prevent the formation of 
a proper torus (i.e. its closing, in the azimuthal direction). In such models, the mass accretion 
rate is only slightly less than the corresponding Bondi rate. 

Subject headings: accretion, accretion discs - black hole physics - galaxies: active 
1. Introduction Narayan & Yi 1994), in which the rate of accre- 



Most galactic nuclei spend a substantial frac- 
tion of their lives in an inactive ("quiescent") 
mode. From theoretical point of view, this inactiv- 
ity is quite surprising because most galaxies, if not 
all, contain a super massive black hole (SMBH) at 
their centers (e.g., Kormendy & Gebhardt 2001) 
and a large amount of gas is available for black 
hole accretion. Thus, one would expect vigorous 
accretion activity resulting in significant emission 
of electromagnetic radiation from all galactic nu- 
clei, but not just from some which are referred to 
as active galactic nuclei (AGN). 



tion within the radius of influence of SMBH can 
be conveniently expressed by the formula derived 
by Bondi (1952). Some other models focus on the 
scenario where the accretion rate itself is much 
smaller than the Bondi value due to rotation, mag- 
netic fields, or both, that can lead to convection 
and mass outflows (e.g., Begelman & Meier 1982; 
Paczyhski & Abramowicz 1982; Narayan & Yi 
1995; Igumenshchev & Abramowicz 1999; Bland- 
ford & Begelman 1999; Stone, Pringle & Begelman 
1999; Quataert & Gruzinov 2000; Machida, Mat- 
sumoto & Mineshige 2001; Hawley & Balbus 2002; 
Proga & Begelman 2003b; Krumholz et al. 2005). 



The modeling of the inactive mode is typi- 
cally based on the assumption of radiatively inef- 
ficient accretion (Ichimaru 1977; Rees et al. 1982; 



In this paper, we focus on exploring effects of 
gas rotation on the black hole accretion hydrody- 
namics, using numerical simulations of gas with 
simplified microphysics . We assume that gas 
accreting onto SMBH is inviscid and its specific 
angular momentum ranges from zero to some fi- 
nite value. Such a range of specific angular mo- 
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mentum is possible for instance in Sgr A* where 
there are many massive stars orbiting the cen- 
tral SMBH (Genzel et al. 2003; Schodel et al. 
2003). As argued by Loeb (2004), winds from 
these stars might be a source of very low angu- 
lar momentum gas. Our study is also relevant 
to other astrophysical situations, e.g., the collap- 
sar model for long duration gamma ray bursts 
(GRBs) where a rotating envelope of an evolved 
massive star collapses onto a central compact ob- 
ject (Woosley 1993; Paczyhski 1998; MacFadyen 
& Woosley 1999; Proga et al. 2003). 

Proga & Begelman (2003a; hereafter PB03) 
studied the axisymmetric hydrodynamical model 
of the slowly rotating gas. Here we generalize this 
model to account for the non-axisymmetric effects. 
This is of a particular interest because the stellar 
winds at the Galactic Center are likely to feed the 
central black hole in a non-axisymmetric way (e.g., 
King et al. 2005; Volonteri et al. 2007). However, 
non-axisymmetric effects can be important even 
for the axisymmetric initial and outer boundary 
conditions because of hydrodynamical instabili- 
ties. PB03 showed that a pressure/rotation sup- 
ported torus forms around the black hole and the 
accretion rate is smaller than the Bondi rate. In 
the inviscid case, accretion is possible only through 
the polar funnels where gas does not rotate or ro- 
tates very slowly. Our goal is to check whether 
this result holds in 3-D and to what extend when 
for example, a non-uniform initial distribution of 
the specific angular momentum is assumed, when 
there is a gas with zero angular momentum at the 
equator. One can expect that if the amount of 
the non-rotating material near the equator is large 
enough, the torus may not form and the accretion 
rate would not be much smaller than the Bondi 
value. One of our main questions is whether the 
accretion can proceed through the equatorial plane 
or it can proceed only through the polar funnels 
as in an axisymmetric case. 

The content of the article is the following. In 
Section [2j we describe the method used in our cal- 
culations. In Section [3] we present the simula- 
tions results. First, we discuss the test runs per- 
formed with the new version of ZEUS-MP code 
(Sec. 3.1l, then, we describe the results for the 
initially axisymmetric 3-D case, which is our 'ref- 
erence' model (Sec. 3.2), and finally we present 



tial conditions (Sec. 3.3 1. We discuss our results 
in Section |U 

2. Method 

In our calculations we use the 3-D code ZEUS- 
MP (Stone & Norman 1992; Hayes & Norman 
2003). The code solves the equations of hydro- 
dynamics: 



dp 



pVv = 



p-^(-)+PVv = 
at p 
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where p is the gas density, e is the internal en- 
ergy density, P = (7 — l)e is the gas pressure, and 
v is the velocity of the flow. We adopt the ra- 
tio of specific heats to be 7 = 5/3. We modified 
the ZEUS-MP code to use the pseudo-Newtonian 
gravitational potential (Paczyhski & Wiita 1980): 



$(r) = - 
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where Rg — 2GM/c 2 is the Schwarzschild radius. 
We start the simulation with the spherical Bondi 
accretion solution (PB03), derived iteratively for 
the density, p(r), energy, e(r), and radial velocity 
v T (r) distributions. We express the accretion rate 
resulting from simulation in the units of the Bondi 
accretion rate: 



M F 



, G 2 M 2 

A47T = poo 



(•5) 



where A 0.29 (see e.g. PB03 for the exact ex- 
pression for A). 

The initial velocity vg is set zero everywhere, 
whereas the velocity v$ is initially non-zero only 
in a fixed, quasi-conical zone, with a limited ra- 
dial size. This zone of initial rotation is formally 
defined as follows: 
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the results for models with non-axisymmetric inl- 
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where l — 0.1 is a dimensionless model param- 
eter. The other model parameters are: Cqo, p^, 
A/bh 1 fin, and r out . The Bondi radius is equal 
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to i?B = GM/c^, and in our simulation it is 
equal to 1000 i?s- Our computational zone ex- 
tends from 1.5xlO~ 3 i?B to 1.2i? B - When present- 
ing the results, we use the units of Rb, Coo, 
and time t — t/t OT \,(RB)- We perform our simula- 
tions for several values of the parameter A<pQ. We 
use the spherical coordinate system, RTP, and the 
boundary conditions in the r, 9 and 4> directions 
are outflow, reflection and periodic, respectively. 
The resolution in r-direction was 140 zones, with 
dri + i/dri = 1.05, in ^-direction it was 96 or 100 
zones, and in (^-direction we had 10, 32 or 60 zones, 
with d8j +1 /d9j = d0k+i/<#k = 1-0. We note here 
that because of a moderate resolution in the </> di- 
rection, the simulations are able to capture only 
the lowest few non-axisymmetric modes. 

3. Results 
3.1. Test runs 

Initially, we performed test runs, to check 
whether the spherical or axial symmetry is con- 
served wherever it should be conserved, and 
whether the 2-D results of the axisymmetric calcu- 
lations are reproduced in 3-D. Having calculated 
both the models of spherical Bondi accretion (s- 
models), as well as the axisymmetric accretion 
with low angular momentum (1-models) , we found 
that: 

(a) the results from the 2-D models are reproduced 
in the r — 6 slices of the 3-D models; in particular, 
the accretion rate M- ln and the flow pattern are 
the same; 

(b) the non-radial velocity components in the s- 
models, which analytically should be vg = iu = 0, 
can locally have non zero values, but both are 
orders of magnitude smaller than the local radial 
velocity and sound speed (that ratio is of the order 
of 10- 10 -10~ 13 ); 

(c) the symmetry with respect to both the equato- 
rial plane and rotation axis is conserved (5x/x < 
10~ 7 , where x denotes density, energy, or velocity 
components) in all models at early stages of the 
evolution, i.e. t < 0.033. 

The s-models conserved their spherical symme- 
try in both 2-D and 3-D, throughout our runs. For 
1-models, as the system evolved, asymmetries be- 
gan to arise with respect to the equatorial plane, 
in both 2-D and 3-D models. We note that this 
was already the case in PB03 simulations. The 



axisymmetry in the present 'reference' 3-D model 
(i.e. the model with axisymmetric initial condi- 
tions) is conserved initially and at intermediate 
times of the system evolution (however see the re- 
sults below). 

3.2. Time evolution of the initially ax- 
isymmetric flow 

The computations of an axisymmetric, 2.5-D 
model of an accretion flow with low angular mo- 
mentum were presented in PB03. We have re- 
calculated this model in 3-D for the purpose of 
the present work. The simulations start from a 
spherically symmetric gas cloud around a black 
hole, with density and velocity distributions de- 
rived from the Bondi solution. The matter located 
far from the black hole possesses a specific an- 
gular momentum that exceeds the critical value, 
Z cl it = 2i?sc, at the equatorial plane and is de- 
creasing towards the polar regions (cf. Eq. 

The time evolution of the system with the ini- 
tial conditions described above proceeds as fol- 
lows. After a transient episode of purely radial in- 
fall, when the rotating material reaches the vicin- 
ity of the black hole, a thick torus forms in the 
equatorial region. The gas settled in this region is 
supported against gravity by the gas pressure and 
rotation, and the rate of accretion on the black 
hole, Mi n , decreased down to about 30% of the 
Bondi accretion rate (in PB03 the exact value of 
M; n was found to depend on the details of the an- 
gular momentum distribution). This is because 
the material accretes only through the polar fun- 
nels, while the torus is made of material that can- 
not accrete. (There is no transport due to viscos- 
ity; however the non-axisymmetric shocks can re- 
sult in some transport of the angular momentum) . 
The gas which approached the centrifugal barrier 
at the equator, could cither outflow radially, or try 
to turn towards one of the poles and accrete. Con- 
sequently, meridional circulation movements are 
observed in the flow. 

The accretion rate onto BH varies in time. Due 
to meridional circulations, the flow is not symmet- 
ric with respect to the equatorial plane; however, 
the time averaged properties (e.g. accretion rate) 
and the shape of the torus did settle down to a 
steady state, as shown in PB03. In the equator 
the flow is subsonic down to very small radii, while 
at the poles, at some distance from the center, the 
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radial velocity exceeds the local speed of sound. 

Here, we repeated the calculations of PB03 in 
3-D and ran new simulations on a relatively long 
time scale. Most of our simulation runs lasted up 
to t = 0.36. For comparison, the sound crossing 
time at the Rb is about i sound (-RB) = 0.125, while 
at the radii corresponding to the dense tori formed 



in the innermost part of the flow (see Section 3.3 1, 
the sound crossing time is * sound (i?torus) ~ 10~ a . 

Figure [l] shows the sonic surface, i.e. the 
isosurface where the Mach number, M to t — 
e + v?) j<? s is a unity, at the t 



0.1. 



For the zero-vorticity flow, the surface shape can 
be derived analytically, as it passes orthogonally 
through the velocity equipotential surfaces (Pa- 
paloizou & Szuszkiewicz 1994). In general, the 
shape of the surface can be more complex. As one 
can see in the Figure, the sound waves are propa- 
gating outwards in the flow. The sonic surface for 
initially axisymmetric model is shown in the left 
panel, and for comparison in the right panel we 
also show the non-axisymmetric case (model A, to 
be described in Sec. 
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We find that both qualitatively and quantita- 
tively the initial results in 3-D axisymmetric case 
are the same as in 2.5-D. The main new feature 
of the torus in 3-D, which could not be studied 
in the PB03 simulation, is that in the very late 
stages of the evolution (t >^ 0.3) the axisymme- 
try breaks and the inner torus becomes tilted with 
respect to the equator, and in the end it starts 
precessing. This departure from the symmetry at 
later times is caused by the equatorial outflow and 
meridional circulations in the torus become sup- 
pressed, while the material tries to get through 
towards the black hole and accrete along one of 
the poles. At the same time, for some regions (e.g., 
4> ~ 0°) the accretion occurs through the northern 
pole, at the opposite side ((f) ~ 180°) and the flow 
chooses rather the southern pole. Consequently, 
a torque is induced and the rotation axis of the 
torus changes in time. 

The precession of the torus is illustrated in 
Figure [2] in terms of the total angular momen- 
tum, itot, which is changing in time. This fig- 
ure shows the motion of the innermost part of the 
flow, i.e. the torus, defined by the density thresh- 
old p m i n = 500. Initially, the dominant component 
is L z , while L x and L y are close to zero (but fluc- 



tuating). It means that the torus rotates basically 
around the z axis. After t ~ 0.21, L x and L y are 
non-zero, and the rotation axis tilts towards the 
4th quarter in the x — y plane. At t ~ 0.3, L y 
starts decreasing while L x still decreases and L z is 
almost constant. It means that the torus is now 
precessing, i.e. the rotation axis moves counter- 
clockwise with respect to an observer along the 
+z axis. Figure [2] also shows for comparison, the 
results for the torus precession in case of non- 
axisymmetric initial conditions (model A). This 
model will be discussed below in more detail (Sec. 

The precession is also shown in Figures [14] and 
15] (see e.g. Fragile et al. 2007). The tilt, defined 
as an angle between the angular momentum vector 
of the gas and the z axis, is initially equal to zero 
for all radii. Later during the simulation the tilt 
rises strongly in the inner parts of the flow. The 
tilt is equal to: 



(3(r,t) 



(t) 



(7) 



In the Figure 14 we plot the tilt angle, (3(r,t), 
as a function of radius, for several time snapshots, 
starting from t = 0.18. The differential form of 
the torus precession is also visible in Figure |15| in 
which we plot the cumulative twist angle, 7(r, t), 
as a function of time. The twist angle is defined 
as a cumulative angle by which the angular mo- 
mentum vector revolves in the x — y plane, by the 
time t: 
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Before the disk was tilted, it did not precess, and 
by definition the twist was zero. Therefore the 
results in the Figure are also plotted from time 
t = 0.18. The solid line is the twist averaged 
over the radius for the inner part of the flow 
(1.5 x 10~ 3 i? B < r < 5 x 10~ 2 i? B ), i-e. the torus, 
while the dashed line shows the twist averaged for 
the whole range of radii. Clearly, the innermost 
torus precesses much stronger than the rest of the 
flow, and the maximum twist at the end of the 
simulation was 7 sa 120°. 

To check whether the tilt and precession of the 
inner torus found in this simulation is a physical 
or rather numerical effect, we performed several 
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further test simulations. First, we reversed the 
direction of the flow in the initial condition, i.e. 
we changed the sign of the azimuthal velocity v^. 
In this simulation, we also found that the torus 
tilts and starts precessing, at the same time (t > 

0. 18). However, the tilt and precession are in the 
opposite directions as measured with respect to 
the grid, i.e. the initial tilt was (3 — 180° and 
decreased to about (3 = 150°, while the twist was 
negative. 

Second, we checked how transient the effect 
of precession is. Due to technical limitations, 
we were not able to run all the simulations for 
very long time, but we completed one run up to 
t = 0.8, for the model R^- In this model we 
found that the tilt angle increased from ~ 30° to 
~ 40°, and the tilt spreaded to larger radii, so 
that not only the innermost parts of the flow pro- 
cessed. However, at this very late phase, some of 
the ring-like structures of the largest density were 
broken into two separate parts, each of them of a 
"C" shape. Therefore the precession may not be 
a long-term effect. We plan to investigate this in 
near future. 

Third, we checked for the importance of the 
adopted physics of the model, in particular, the ra- 
tio of the sound speed to the free fall velocity. We 
calculated two models with much smaller Bondi 
radius, Rq — 100i?g and R^ = 300i?s, which cor- 
respond to a much larger sound speed at infinity: 
respectively, 3.15 and 1.82 times larger than in all 
the other simulations. These models are denoted 
as Mioo and Af3oo in Table [l] The models could 
be run for much longer time in terms of i rb(-RB), 

1. e. for t — 11 and t — 2, respectively. How- 
ever, in these models we did not find any signa- 
tures of tilt or precession. This is because here the 
Mach numbers are never large: they are at most 
M = 3.3 — 3.5 at the inner radius, whereas for the 
precessing torus the Mach numbers reached the 
values as large as 4.7 - 6.0. The small Mach num- 
bers in models Mioo and M300 make the shocks 
smaller; hence they do not amplify the asymme- 
tries growing out of initial perturbations. 

Finally, we tested the role of the artificial vis- 
cosity, which might help to spread the shocks and 
avoid precession if it was a numerical artifact. The 
artificial viscosity was parametrized with the stan- 
dard Neumann-Richtmeyer artificial viscosity co- 
efficient qcon =2.0 However, the results in this 



simulation were very similar to the original sim- 
ulation. Specifically, for time t — 0.324 the max- 
imum tilt angle was 26°, while in the former case 
and 28° in the latter. 

From the above tests, we conclude that the 
precession in our initially axisymmetric model is 
rather a physical than numerical effect, and is con- 
nected with relatively large supersonic speeds of 
the flow achieved in our model. 

3.3. Time evolution for the non-axisymmetric 
initial conditions 

Now we investigate how the non-axisymmetric 
initial distribution of the specific angular momen- 
tum affects the evolution of the flow. In particu- 
lar, we check whether the rotationally supported 
torus forms and if a steady state can be achieved 
(with or without torus precession). We start our 
simulation with the non-zero specific angular mo- 
mentum enclosed in a conical region of a width 
A(j)Q (see Eq. [6]). A naive prediction could be 
that the rotating gas will reach the innermost re- 
gions, spiral in, and after a few orbital cycles the 
material with large angular momentum would be 
mixed with the non-rotating gas. Therefore a ro- 
tationally supported torus would form, regardless 
of A <po- The only dependence on this parameter 
would be the moment when such a torus forms. 
However, as we show below, the numerical simu- 
lations lead us to a different result: depending on 
A(f>Q the rotating gas may not form a torus at all. 

We performed the runs for non-axisymmetric 
initial conditions for a range of A</> . The models 
are summarized in Table[T] The non-axisymmetric 
models are labeled with the letters A-E, while 
the reference model is labeled as R. As we men- 
tioned in Sec. [2] and as the Table shows, we tested 
the models with smaller (32 zones) and larger (60 
zones) resolution in the ^-direction. We checked, 
that the time averaged results for the accretion 
rate only very weakly depend on the resolution, 
however the amplitude of time variability of M 
increases with resolution. 

Below, we present these results for the largest 
adopted value of A</>o = 330° (model ^4). This pa- 
rameter translates into a small non-axisymmetric 
perturbation in the initial conditions, i.e. small 
content of non-rotating material. 

In Figures [3] [4] [5j [6] and [7] we show the color 
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coded maps of the central region, as well as the 
zoomed-out. The maps show the density distri- 
bution and velocity field, as well as the specific 
angular momentum, plotted for several snapshots 
during the evolution: i =0, 0.018, 0.09, 0.16, 0.23 
and 0.29. Note that the top-right panel in Figure 
[3] as well as all the top panels in Figure [6j i.e. 
for l spec (t = 0), are plotted on the scale 20 times 
larges compared to other panels, because in the 
inner region initially we assumed l spC c = 0. 

Figure [3] shows the density and specific angu- 
lar momentum distribution in the central region 
(up to 0.02 Rq), plotted for the equatorial plane. 
The orientation of the plots in the x — y plane is 
standard, i.e. the x > semi-axis corresponds to 
= 0. The density distribution, initially spher- 
ical at t = (top left panel) changes in time, as 
the material which carries specific angular momen- 
tum approaches the center. At t =0.018, the gas 
is rotating around the z-axis (arrows over plot- 
ted on the density maps denote the direction of 
the velocity vectors with components v v and v^,). 
The material is distributed axisymmetrically, and 
the specific angular momentum is rather large (i.e. 
0.8 < /spcc/^crit < 1-2) for most of the 4> directions. 
A clump of gas with relatively smaller l spQC , which 
can be seen on the second right panel, is a remain- 
ing of the gas with / spcc = initially present at the 
equator for <fi = (—15°, 15°). This clump is track- 
ing the archimedean spiral (see also e.g. Lemaster 
et al. 2007), and such a trajectory appears to be 
due to the pressure gradient force and rotation 
(similarly to a cyclone) . 

At t = 0.09, the material with small Z spcc is al- 
ready mixed with the gas of high Z spcc . The mate- 
rial rotates very fast in the equatorial plane, and 
a circular pattern of slightly larger and smaller 
specific angular momentum regions, visible on the 
third plot, form due to mixing. These motions 
are suppressed at t =0.16. For some specific di- 
rections, namely (f> ~ 90° and c/> ~ 270°, the angu- 
lar momentum near the inner radius becomes very 
small (as indicated by the green spots in Fig. 
and the density in these regions drops, indicating 
that the material falls radially to the black hole. 

The solid line in the angular momentum maps 
mark the contour at which Z sp cc/fcrit = 1-0. We 
note that at t' = 0.23 and t — 0.29, in the equa- 
torial plane a substantial fraction of material has 
^spoc < ^crit- This material is located at an ex- 



tended region approximately along the diagonal of 
the plane (i.e. <f> ~ 45° and 4> ~ 225°, as marked 
by the orange shade in Fig. |3j. Another region of 
even smaller angular momentum, Z spcc <C Z c rit> has 
a very limited radial extension close to the center, 
and is slightly elongated along the other diagonal 
(i.e. 4> ~ 135° and <j> ~ 315°, as marked by the 
green shade in Fig. |3j. For the same 0, but at 
larger distances, the specific angular momentum 
in the equatorial plane is very large, and exceeds 

Zcrit ■ 

In Figure [4] we show the distribution of l spcc on 
the larger scale. (Note also that the color scale is 
different than in Fig. [3] The red and green shades 
now mark the regions with Z S pec/Z C rit > 1.0). As 
the Figure shows, more material with l spQC < Z cr it 
mixes in at late stages of evolution, t = 0.09 sec. 

From the Figs. [3] and [i] we conclude that a 
torus starts forming in the innermost regions in 
the equatorial plane at t =0.018, and it is ro- 
tationally supported for all the <f> directions, be- 
cause the low Z spoc material is quickly mixed in. 
At t =0.09, the distribution of density and spe- 
cific angular momentum is almost axisymmetric, 
but only in the inner region. However as the zoom- 
out maps show, the gas with Z spcc < Z cr it is mix- 
ing in and is distributed asymmetrically The in- 
nermost symmetric configuration lasts until about 
t =0.16, and after this time the torus position 
and shape changes. To investigate what really 
happens, we need to look at the flow from a dif- 
ferent perspective. Therefore, in Figures [U|6] and 
[7] we show the slices perpendicular to the equa- 
torial plane. The maps show density and velocity 
field, as well as the specific angular momentum, as 
seen from cj> = 0°, 90°, 180° and 270°, at the same 
times as in Figs |3] and [4] 

As shown in the Figure [E] at i =0.018, the 
material is indeed accumulated near the equator, 
while at the poles the density is much lower. The 
axial symmetry is not perfect, since at = 180° 
the torus is geometrically thicker than at other <j> 
directions, which corresponds to the location of 
the 'clump' with smaller specific angular momen- 
tum (c.f., Fig. [3]). As indicated by the arrows 
(velocity vectors with v r and vg components), an 
equatorial outflow occurs at most of the <$> direc- 
tions, however at <f> = 180° this outflow is weaker, 
because of slower rotation. 
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At t = 0.09, the torus is relatively thin and lo- 
cated at the equator in every direction, while at 
the poles the density is very low. The gas accretes 
onto the center mostly through the poles, while at 
the equator the flow pattern is complex (circula- 
tions, outflows). The flow complexity is reduced at 
t = 0.16, because the gas that flows into the equa- 
torial region has less angular momentum and can 
directly accrete onto the black hole, finding its way 
along one of the poles. At <f> = 0° and <j> — 270° 
the gas turns towards the northern pole (i.e. +z), 
while at </> = 180° and cj) = 90° the flow turns to- 
wards the southern pole (i.e. — z). A line which 
would mark the regions of the maximum density, 
is now tilted with respect to the equator by an an- 
gle of ~ 20°. This means that the torus which at 
t — 0.09, was almost symmetric with respect to 
the equatorial plane, is tilted after t — 0.16. 

In Figure [6] we show the maps of specific angu- 
lar momentum, also perpendicularly to the equa- 
torial plane and for the same <$> directions as in 
Fig. [5] As the Figure shows, the specific angular 
momentum distribution is symmetric with respect 
to the equatorial plane at t = 0.09, and nearly 
symmetric at t = 0.16. The axial symmetry is 
not perfect, however the regions of large angular 
momentum (i Bpec > i C rit) appear in every slice. At 
t > 0.16, the flow is no longer symmetric, and the 
regions with very large specific angular momentum 
appear either below the equator (at <f> = 0° and 
(f> = 90°) or above it (at 4> = 180° and 4> = 270°). 
These regions correspond to a flatter torus, i.e. 
relatively thin on one side, while the gas which 
is not rotating fast makes the configuration geo- 
metrically thicker. We notice, that the lack of the 
top-bottom symmetry in density and l spC c maps at 
this phase of system evolution is a consequence of 
the earlier non-axisymmetry. This asymmetry was 
introduced to the velocity field in the initial condi- 
tions at large radii, and subsequently propagated 
to the inner radii and affected the density distri- 
bution there, as soon as the rotating gas reached 
there. 

In Figure [7] we show the distribution of the spe- 
cific angular momentum in the zoom out. The 
flow is rotationally supported in the outer regions. 
The Z spoc distribution in the flow is asymmetric at 
large scales (there is significantly more material 



the case only in the inner region, i.e. the torus is 
tilted, while outer regions, even at late times, are 
rather symmetric with respect to the equatorial 
plane. 

To visualize the 3-D configuration better, in 
Figure [8] we show the density isosurfaces in 3-D. 
The plots show 3 arbitrarily chosen contours of the 
constant density: 2000, 1250 and 500 Poo, which 
correspond to the gas densities very close to the 
inner radius, i.e. inside ~ 0.004i?e- The maps are 
plotted for 5 different time snapshots: t =0.018, 
0.09, 0.16, 0.23 and 0.29. We do not show the 
density distribution at t=0, because it is purely 
spherical. The orientation of the figures is almost 
edge-on, i.e. the z-axis is the rotation axis of the 
system, and x — y plane is the equatorial plane. 

These density contours may be regarded as the 
shapes of the torus (however one should keep in 
mind that the material with smaller/larger den- 
sity is also present there). Therefore, as the Fig- 
ure shows, the torus is closed already at t = 0.018. 
After t — 0.16, the configuration becomes tilted 
with respect to the equatorial plane, and at t = 
0.29 this tilt is the largest. After t = 0.29 the 
torus starts precessing, and the precession period 
is very long (by the end of the run, the ring pro- 
cessed around the z-axis by less than 90°). The 
torus precession was also shown in the right panel 
of the Figure [2] (Sec. 3.2 1. The changing values 
of L x , L y and L z , show that at t ~ 0.1 the rota- 
tion axis, which is defined by the direction of L, 
starts to tilt towards the first quarter of the x — y 
plane, and after t ~ 0.29 the rotation axis moves 
counter-clockwise . 

In Figure [9] we show the isosurfaces of the spe- 
cific angular momentum in 3-D. The contours are 
for i/Z cr it = 1.3, 2.15 and 3.0, and correspond to 
the zoomed-out regions shown in Fig. [7] (the ra- 
dial extension of about 0.6 Rb)- As the Figure 
shows, the material with various angular momen- 
tum is mixing in the innermost region at t — 0.018 
while at the outer parts initially the momentum is 
not mixed. At later times, the outer parts of the 
flow contain more and more mixed angular mo- 
mentum layers, while in the innermost region the 
gas rotates slower than at the beginning, and the 
angular momentum distribution is smoother. 



with large Z spec for the directions of 



0° and 



270° than for 4> = 90° and 180°). However, this is 



In Figure 10 we show the maps of entropy, S, 
radial to azimuthal velocity ratio, v r /v^,, angular 
velocity, f2, and velocity divergence, div v, as cal- 
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culated close to the end of this simulation, at time 
t =0.29. The maps show the inner region in 
the equatorial plane. Close to the center, the an- 
gular velocity is the largest along the diagonal of 
the plane (<f> ~ 45° and <fi ~ 225°), which cor- 
responds to the cross-section line along which the 
torus crosses through the equatorial plane (cf. Fig. 
[8]). This line also corresponds to the largest en- 
tropy, as well as positive velocity divergence, while 
at the other diagonal there is smaller entropy and 
negative velocity divergence. This means that the 
fluid is compressible (and supersonic). 

Our analysis of the 3-D results shows the flow 
has a largest negative divergence close to the poles 
while the largest positive divergence is somewhat 
above and below the equator, i.e. on the surface 
of the torus, as well as at its cusp. From the poles, 
the gas flows radially onto the center with large su- 
personic velocities, i.e. M T = yf /c% ^> 1. Close 
to the equator, the flow is captured in the torus 
and the radial velocities are smaller, so M r < 1. 
Still, the gas rotates very fast, and the total Mach 
number is large, M to t 3> 1> because of the contri- 
bution from the azimuthal velocity. 

The entropy in the flow should be constant 
along the streamlines, however may vary from one 
streamline to the other. The entropy gradient cor- 
responds to a non-zero vorticity in the flow. The 
direction of vorticity arrows (w = rot x v) in the 
bottom-right panel indicate that the torus is ro- 
tating counter-clockwise. The direction of radial 
velocity arrows in the bottom-left panel confirms 
that the radial inflow occurs from the directions 
of the smallest angular velocity. 

The above considerations led us to the conclu- 
sion that qualitatively, the basic pattern of the 
torus evolution is uniform. It does only weakly 
depend on whether we assume the axisymmetric 
initial conditions, or if the initial distribution of 
angular momentum was perturbed (provided that 
this perturbation was small enough to allow for the 
torus formation). The axisymmetry was explicitly 
assumed by PB03 throughout their 2-D simula- 
tions, as well as in our 3-D reference models in the 
initial conditions. The similarity of the torus be- 
havior that we find here means that in both 3-D 
models R and in the A models, the rotationally 
supported tori form after time t = 0.018, than 
they exhibit strong equatorial outflows (which sta- 
bilize them) , then they become tilted with respect 



to the equatorial plane due to the asymmetric po- 
lar accretion, and finally start precessing. What 
differs in the models, is the moment when the out- 
flow stops and when the torus becomes tilted, as 
well as the tilt angle (it is about twice as large in 
model A than in model R at the end of our simu- 
lations; see also Fig. j2|. These two features (i.e. 
tilt and precession) cannot be studied in the 2-D 
simulations, but are detected in the 3-D models. 

However, we note that the similarity between 
the R and A models is less pronounced when con- 
sidering the details, e.g., of the torus shape. The 
asymmetric torus can be thicker (warped) , or thin- 
ner, depending on (f>. To investigate further the 
difference between R and A models and the role of 
non-axisymmetry, we calculated the angular mo- 
mentum at the inner boundary, as a function of 
the angles 9 and 4>. 

In the axisymmetric model R, the specific an- 
gular momentum at r m depends only on 9 by defi- 
nition. After the torus is formed, the maximum 
value at the equator is Z S pecAcrit(?* = Hn,^ = 
90°) « 0.85, while the value averaged over the 
angle 9 is l(r — r; n ) w 0.55. For the non- 
axisymmetric case (in particular, model A), the 
results depend also on the angle 4>. Therefore both 
the equatorial angular momentum, l(r = r m ,9 = 
90°), and the 6>-averaged, l(r — r in ), are scattered 
and do not have to match with the axisymmetric 
solution. In Figure [TT] the results for the axisym- 
metric model, are denoted by single points, while 
the non-axisymmetric solutions are represented by 
the horizontal lines to show the scatter in <j). 

The level of non-axisymmetry is represented as 
the spread between the maximum and minimum 
values of Z spG c at a given time. For the model 
A, in the beginning of the torus evolution, i.c 
0.018 < t < 0.09, the equatorial and averaged 
values of Z spoc match well with the axisymmetric 
solutions and are only slightly scattered with cf>. 
The equatorial value is much larger than the 9- 
averaged, which means that the torus is located 
in the equatorial plane where ^ spcc is the largest. 
As the evolution proceeds, t > 0.09, the scat- 
ter with <j) increases and the range of l apcc does 
not match the axially symmetric solution. How- 
ever, the equatorial angular momentum is still 
much larger than the average. This means that 
the torus is located at the equatorial plane, but 
it is now asymmetric. After t > 0.16, the situa- 
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tion changes, and l spec at the equator has a very 
large scatter, being either smaller or larger than 
the 0-averaged value (i.e. the solid and dashed 
lines overlap in the Figure). The averaged angu- 
lar momentum has also some scatter, but much 
smaller than the equatorial one. This indicates 
that the torus is not symmetric, and is tilted with 
respect to the equatorial plane. 

We performed similar analysis of other non- 
axisymmetric models, B-E (cf. Table [lj. For 
models B and C (A0 O = 240° and A</> = 120°), 
at the beginning of the simulation 0.018 < t < 
0.07, the equatorial value of Z spcc is always larger 
than the averaged, and the scatter with <fi is 
smaller in model B than that in model C. This 
indicates that a rotationally supported torus is 
present in the equatorial plane. At later times, 
the equatorial Z spcc becomes equal or smaller than 
average, and the scatter in both models B and C 
is quite substantial (larger than in model A) . This 
implies that a torus, which possibly tried to form 
at the early phase, is broken (i.e. not completely 
closed), as well as tilted from the equator. An 
example of such a 'broken torus' configuration is 
shown in Figure |16| 

For models D and E (A<j) = 30° and A0 O = 
60°), the scatter in both 0-averaged and equato- 
rial l spcc is very large and does not decrease with 
time (up to t ~ 0.018). The averaged l spcc can 
be larger or smaller than the axisymmetric one, 
while the equatorial Z spcc in these models is al- 
ways smaller than that in the axisymmetric case, 
and locally (i.e. for some <f> angles) can be smaller 
than the 6- averaged. We conclude that in these 
models the torus does not form, the solution is 
not axisymmetric, and the gas with very small an- 
gular momentum can accrete onto the black hole 
through the equator. 

We note that in this sense the properties of 
these models are similar to those of the model A 
in later times. However, when comparing the den- 
sity distributions, in the models D and E the gas 
is always distributed much more uniformly, i.e. it 
does not concentrate neither close to the equator 
nor to any specific plane. In model A. the gas 
density near the poles is always orders of magni- 
tude lower than that at the equatorial plane, or 
the plane tilted to the equator by a small angle. 

This is not the case for models D and E, for 



which the torus does not form. For models B and 
C, the gas is concentrating near the equator, but 
only for some range of (/(-directions, i.e. the torus 
is not closed. 



Figure 12 shows the time evolution of the ac- 
cretion rate through the inner boundary, Mj n for 
the non-axisymmetric models. Before the rotating 
material approaches the black hole (t < 0.018), 
Mi n is equal to the Bondi accretion rate for all 
models. Once the gas starts rotating also in the 
innermost parts, the accretion rate drops reaching 
about 25% - 40% of the Bondi rate for models A, 
B, C and R. This is because much more mate- 
rial gets captured in the rotating torus, and does 
not fall radially into the black hole. The model A 
gives the lowest accretion rate with rather small 
and regular variability pattern, very close to that 
obtained in the reference model R. 

Figure [13] presents the evolution of the angu- 
lar momentum flux through the inner boundary: 
Lhi — J hpecPVrds, in units of the critical angu- 
lar momentum Z cr it and renormalized by the value 
of the Bondi accretion rate. At the beginning of 
the simulation, L m — 0, since there is no rota- 
tion in the vicinity of the black hole. Once the 
rotating matter reaches the inner boundary, the 
angular momentum starts accreting to the center. 
When the torus starts forming, the fast rotation 
near the center leads to a fast rise in L in . However, 
after several orbital cycles the outflow begins and 
the net radial velocity drops, as well as drops the 
density near the polar regions, therefore the L 1TL is 
rather small during the torus evolution. 

Note that the quantity plotted in the Figure [13] 
is not a flux of specific angular momentum, but the 
total one. This corresponds to the amount of an- 
gular momentum which may be transferred to the 
black hole and used to spin it up. However, as the 
Figure shows, the total angular momentum which 
the black hole could gain during our simulation, is 
extremely low: a = (cJ)/GM 2 « 4 x 10~ 6 - 10~ 5 , 
where J = L m At. 

Figure [13] shows two trends in the magnitude 
of L ln . For small A^o (models E and D), the ro- 
tation at inner boundary is very small, but the 
density is high both in the equator and in the po- 
lar regions (still close to the spherical accretion). 
Therefore, in model D, with faster rotation, Lj n is 
larger. For large A</>o (models A, B and C), the 
material accumulates rather close to the equator, 
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at least for some <j> angles, while at the poles the 
density is small. Therefore regardless of the fast 
rotation, these models give systematically smaller 
Lin than those in the models E and D. Also, this 
is why the model B gives smaller L ln than that in 
the model C. Model A is the only one in which the 
torus is closed (the gas rotates fast at the equator 
at every <j> angle). This leads to a larger L in than 
that in model B, which is again less affected by 
the density distribution. 

In the models A, B and C, the accretion rate 
and Lin are variable. For model A, a characteris- 
tic wave pattern can be seen in the specific angu- 
lar momentum distribution in the equatorial plane 
(see Fig. |i]at t ~ 0.09). This behavior is reflected 
in the variable Mi n and L- ln , due to the variable 
radial velocity and nearly constant density at the 
inner boundary. At time t ~ 0.16, when the ac- 
cretion rate stops varying rapidly, the correspond- 
ing 'waves' in specific angular momentum map are 
smoothed out. The outflow of the gas at this time 
is suppressed, and the gas accretes onto the cen- 
ter through the poles. As a consequence, the mass 
flux and angular momentum flux through the in- 
ner boundary slightly increase, and the curves are 
smoother. The subsequent drop of both of these 
quantities at time t = 0.25 is caused by the den- 
sity decrease at the inner edge, when the torus 
is tilted with respect to the equatorial plane. The 
density at inner radius increases again at t — 0.29 
and the torus starts precessing. 

The behavior of the flow in models B and C is 
more chaotic. The accretion rate and L m vary in 
time until the end of our simulations. Also, no pre- 
cession was detected, because in principle it would 
be hard to determine the boundaries of the struc- 
ture which might be precessing. We cannot use a 
density threshold to define a torus, because such 
a torus is not a closed ring in these two models. 

4. Discussion 

This paper presented the fourth phase of our 
study of slightly rotating accretion flows onto 
black holes (see PB03, Proga & Begelman 2003b, 
and Proga 2005 for the first, second and third 
phase). Here we followed PB03, but we consid- 
ered 3-D not 2-D axisymmetric effects. As in 
PB03 we made a few simplifications. For exam- 
ple, we neglected the gravitational field due to the 



host galaxy, radiative heating and cooling, vis- 
cosity and MHD effects. Perhaps the most im- 
portant simplification we made is neglecting the 
transport of energy and angular momentum out- 
ward as needed to accrete matter with a specific 
angular momentum higher than l cr n- As shown 
Proga & Begelman (2003b) for the axisymmetric 
case, magnetic fields, that can drive the transport, 
can dramatically alter the flow solution. Still our 
HD results provide a useful exploratory study of 
accretion onto black holes as they have revealed 
unexpected properties and complexity of accretion 
flows even with simplified physics. To our best 
knowledge there have not been any MHD simu- 
lations of rotating flows where a closed torus was 
cither assumed or failed to form. Our results show 
that such flows are plausible (see also Loeb 2004) 
and motivate new simulations. In what follows we 
summarize and discuss our results. 

We have performed numerical 3-D hydrody- 
namical simulations of slightly rotating, inviscid 
accretion flows onto a black hole. As in PB03, 
we attempt to mimic the boundary conditions of 
classic Bondi accretion flows with the only modifi- 
cations being the introduction of a small, latitude- 
dependent and also azimuth-dependent angular 
momentum at the outer boundary and a pseudo- 
Newtonian gravitational potential. The adopted 
form of the distribution of l apec the density distri- 
bution at infinity to approach spherical symmetry, 
because the centrifugal force is negligible. 

For the latitude-dependent angular momen- 
tum, 3-D simulations confirm axisymmetric re- 
sults. Namely, the material that has too much 
angular momentum to be accreted forms a thick 
torus near the equator. Therefore the geometry 
of the polar funnel, where material is accreted, 
and the mass accretion rate through it are con- 
strained by the size and shape of the torus but 
by the outer conditions. However, in 3-D the 
torus precesses and is non-axisymmetric even for 
axisymmetric conditions at large radii. For the 
latitude- and azimuth-dependent angular momen- 
tum in the initial conditions, the non-rotating gas 
near the equator can significantly affect the evo- 
lution of the rotating gas. It can prevent closing, 
in the azimuthal direction, of the rotating gas and 
the proper torus does not form. In such cases, 
the mass accretion rate is only slightly less than 
Bondi rates. 
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Simulations with none or a small amount of 
a non-rotating gas near the equator show that a 
torus forms and limits the accretion rate. How- 
ever, a non-rotating gas near the equator can in- 
hibit torus formation and the accretion rate will 
be close to the Bondi rate. Thus, our simulations 
show that in 3-D it may be even more difficult than 
in 2-D to explain the inactive mode of accretion. 
However, if the torus forms then our simulations 
show that the torus will precess. This precession 
may have important consequences in terms of re- 
ducing the mass accretion rate. Namely, a process- 
ing torus may produce a precessing jet/wind that 
will then affect a larger volume of the surrounding 
material than a non-precessing jet/wind. 

Our simulations are in principle relevant to any 
type of a black hole, which accretes gas with some 
small angular momentum, because the results 
should scale with the black hole mass. Note, that 
the parameters in our models were chosen so that 
the ratio of the Bondi radius to the Schwarzschild 
radius was equal to 10 3 and the computational 
domain was between 1.5 i?s and 1.2 Rb- In this 
context, one of the most interesting cases is that 
of the low luminosity active galaxies, and in par- 
ticular the Sgr A*, for which the supermassive 
black hole of M ~ 3.5 x 1O 6 M was identified in 
the center (Ghez et al. 2005). 

The Bondi accretion rate in the Galaxy center 
inferred from the studies of stellar winds was es- 
timated to be between 10~ 6 and 10~ 4 A'/ Q yr _1 , 
while the polarization studies suggest that the 
actual accretion rate is between 10~ 7 and 10~ 5 
M yr _1 from ROSAT observations (Quataert, 
Narayan & Rcid 1999; Baganoff et al. 2003; Bower 
et al. 2005). Although the values overlap, most 
of the authors agree that in Sgr A* the accretion 
rate onto the black hole is well below the Bondi 
value. Also, studies of other quiescent AGN sug- 
gest that some modification of the Bondi accretion 
is needed (e.g. Di Matteo et al. 2000). 

The hydrodynamical studies of black hole ac- 
cretion in Sgr A* which took into account a ran- 
dom distribution of stars were presented in Coker 
& Melia (1997). The analytical estimates of the 
effective angular momentum of the accreting gas 
performed by Moscibrodzka, Das & Czerny (2006) 
were based on the strengths of the stellar winds. 
Recently, Cuadra, Nayakshin & Martins (2006) 
modeled the wind accretion allowing the stars 



to move on the elliptical orbits. We show here 
that the accretion rate of less than 40% of the 
Bondi value is possible also in models, where 
the rotationally supported torus does not close, 
and an asymmetric structure of material being 
almost spherical on one side of the black hole 
may persist as a kind of steady-state (model with 
Ac/)q = 120°). Only for very small input of the 
angular momentum gas, corresponding to the to- 
tal Ac/)q < 60°, the accretion flow is still remains 
almost spherical and the accretion rate does not 
drop much below the Bondi value. 

Perhaps our most intriguing results we found, 
is the instability and precession of the torus. The 
precession occurs for the closed rotationally sup- 
ported torus, which forms for large angular mo- 
mentum input 330° < A0 O < 360°. In other 
words, we found that even a very small asymmetry 
in the angular momentum distribution (not nec- 
essarily in the initial conditions) will lead to the 
torus misplacement from the equatorial plane and 
its precession. This happens after a few tens of 
orbital cycles. For the models where the torus did 
not form due to a very large content of non ro- 
tating gas accreting in purely radial direction, the 
asymmetric condition makes the quasi-steady con- 
figuration very unstable, and rapid fluctuations of 
the global flow pattern occur on the timescales of 
a few dynamical cycles. The clumps of gas with 
large density become misplaced from the equator. 
For example, in the model with A<j>o = 240° the 
misplacement reaches a few tens of degrees and 
changes in the flow configuration are extremely vi- 
olent. 

The problem of the stability of accretion tori 
with respect to the axial perturbations was stud- 
ied in a number of papers. The classical Raylcigh 
condition for the torus stability (sufficient only for 
axisymmetric modes) is that the specific angular 
momentum should not decrease outwards (see e.g. 
Chandrasekhar 1961). The Kelvin- Hclmholtz in- 
stability occurs when two superimposed layers of 
fluid arc in a relative motion. When the velocity 
shear exceeds a critical value, the resulting pres- 
sure gradient (from Bernoulli's law) between the 
peaks and troughs of an interfacial wave overcomes 
the surface tension and gravity, and the mode 
grows exponentially. 

Papaloizou & Pringle (1984) studied the stabil- 
ity of the non-axisymmetric modes of the differen- 
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tially rotating tori. In this first paper, they lim- 
ited their considerations to the homentropic tori 
with constant specific angular momentum, and 
they found, that all such tori arc unstable to the 
low order modes and the instability occurs on a 
dynamical timescale. These instabilities are found 
to be global, i.e. their presence cannot be detected 
via the local analysis nor from the considerations 
of axisymmetric modes. In their second paper (Pa- 
paloizou & Pringle 1985), they considered the tori 
with non-constant specific angular momentum and 
they found that for low azimuthal number m the 
instability is driven by a Kclvin-Hclmholtz mech- 
anism. The modes are stable in disks with angu- 
lar velocity decreasing with radius as oc r~ q , if 
q > \/3. For high m, the modes regain their sonic 
character, and there exist sonic modes which are 
driven by both mechanisms. However, their ana- 
lytical calculations were done in the limit of a Ke- 
plerian disk (q = 3/2). In our models, q is about 
1.8. In the innermost parts, the flow is supersonic 
and compressible. Future work is needed to per- 
form stability analysis for such flow properties. 

From the point of view of numerical, multi- 
dimensional simulations, as it has been recently 
discussed by Foglizzo, Galletti & Ruffert (2005) 
that it is not always obvious to determine whether 
the hydrodynamical instabilities are a physical or 
numerical effect. As shown by our simulations, 
and confirmed by the tests with smaller Mach 
numbers (models M wo and M 300 ), the instabili- 
ties that lead to the torus precession do not de- 
velop in a subsonic flow (see also Moscibrodzka & 
Proga 2008). The acoustic instability develops in 
the 3-D accretion when the Mach number is large, 
and for 7 = 5/3, and it is found to be a physical 
instability in a strongly supersonic flow. We plan 
to verify this result for other values of 7, and per- 
form more detailed resolution tests in the future 
work. 

From the observational point of view, the pre- 
cessing torus might be relevant to the interpre- 
tation of jet emitting sources, provided that the 
jet axis is always perpendicular to the disk sur- 
faces. Recent observations of radio jets show jet 
reorientation. One of possible explanation for this 
reorientation may be the jet precession, as was 
suggested e.g., for the shape of the source 3C294 
(Erlund et al. 2006). Also, the morphology of 
the BAL quasar 1045+352 indicates either a pro- 



cessing jet or an ongoing merger process (Kunert- 
Bajraszewska & Marecki 2007). In addition, the 
morphology of some of the VLBA observed curved 
jets suggests jet precession because of relatively 
short timescales (Lister 2006). 
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Tabic 1: Summary of the models for non- 
axisymmetric accretion hydrodynamics in 3-D 



Model 


A0 O 


Resolution 


^end 


M 








[°] 


[N T X Ng x N^] 


[t rb{R&)] 


[Mb] 


Presence of torus 


Precession 




330 


140x96x32 


0.36 


0.25 


yes 


yes 


-460 


330 


140x96x60 


0.20 


0.25 


yes 


yes 


-B60 


240 


140x96x60 


0.30 


0.28 


not closed 




C32 


120 


140x96x32 


0.36 


0.37 


not closed 




Ceo 


120 


140x96x60 


0.16 


0.35 


not closed 




A32 


60 


140x96x32 


0.18 


0.90 


no 




-E32 


30 


140x96x32 


0.14 


0.98 


no 




-Rio 


360 


140x100x10 


0.32 


0.24 


yes 


no 


^?32 


360 


140x96x32 


0.36 


0.24 


yes 


yes 


-^60 


360 


140x96x60 


0.36 


0.24 


yes 


yes 


M100 


360 


140x96x32 


11 




yes 


no 


M 300 


360 


140x96x32 


2.19 




yes 


no 
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Fig. 1. — The sonic surfaces for the 3-D mod- 
els with axisymmetric initial conditions (left) and 
the non- axis ym met r ic case (right; model A 32 de- 
scribed in Sec. 3.3 1 at time t =0.1. The ex- 
tension of the surface in z direction is about 0.065 
i?B, while in x and y direction the box size is about 
0.03 R B - 



Fig. 2. — The total angular momentum evolu- 
tion in the initially axisymmetric (left) and non- 
axisymmetric (right; model A32 described in Sec. 
3.3) models. The angular momentum is in the 
units of Keplerian angular momentum on the in- 
ner radius (Zk, and the plot shows the magnitude 
of I/tot (solid line) and its components: L x (short 
dashed line), L y (long dashed line) and L x (dot 
dashed line). 
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Fig. 3. — The results for model A32(see Tab. [IJ, 
for 6 time snapshots, from top to bottom: t = 
0.0, 0.018, 0.09, 0.16, 0.23 and 0.29. The maps 
show the central region (only upper right map is 
zoomed-out) in the equatorial plane (6 = 90°): 
density and velocity field (left) and specific angu- 
lar momentum (right). The solid line marks the 
contour of l spcc = Z crit . 
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Fig. 4. — The same as in Fig. [3] but the maps 
show the equatorial plane in a zoom-out (i.e. 20 
times larger scale). Note that the color scales are 
now different than in Fig. [3] 




Fig. 5. — The results for model A32 (see Tab. [TJ, 
plotted in the r — 6 plane, at 4 slices of the 4> 
angle: 0°, 90°, 180° and 270°. The maps show 
density and velocity fields in the central region. 
The corresponding times are, from top to bottom: 
t = 0.0, 0.018, 0.09, 0.16, 0.23 and 0.29. 
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Fig. 6. — The maps of the specific angular mo- 
mentum for model A32 (see Tab. [T]), plotted in 
the r — 9 plane, at 4 slices of the 4> angle: 0°, 
90°, 180° and 270°, and at t = 0.0, 0.018, 0.09, 
0.16, 0.23 and 0.29 (from top to bottom). The 
maps show the central region (only first map is in 
zoomcd-out). The solid line mark the contour of 

^spcc — ^crit ■ 




Fig. 7. — The same as in Fig. [6] but all the maps 
are in zoomed-out. Note that the color scale is 
now different than in Fig. [6] 
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Fig. 8. — Density isosurfaces for five time snap- 
shots: t = 0.018, 0.09, 0.16, 0.23 and 0.29, for 
model A 32 (see Tab. [TJ. The colors mark the sur- 
faces of p = 500, 1250 and 2000 p^. The boxes 
are scaled to the density isocontour p = 500 p^ 
and show the central region of the radius about 
0.004 R B . 



Fig. 9. — Isosurfaces of the specific angular mo- 
mentum in five time snapshots: t = 0.018, 0.09, 
0.16, 0.23 and 0.29, for model A 32 (see Tab. [TJ. 
The colors mark the surfaces of ^ S p OC /^crit = 1-3, 
2.15 and 3.0. The boxes are scaled to the isocon- 
tour of Z S pec/fcrit =1.3 and show the region of the 
radius about 0.6 Rb- 
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Fig. 10.— The results for model A 32 (see Tab. [T| 
plotted in the equatorial plane, at t = 0.29. The 
maps show: angular velocity (upper left), veloc- 
ity divergence (upper right), radial to azimuthal 
velocity ratio (bottom left) and entropy (bottom 
right). The arrows denote the directions of ve- 
locity vectors with v r and components (upper 
two panels), with only v r component (bottom left 
panel) or vorticity vectors with w r and com- 
ponents (bottom right panel). 
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Fig. 11. — The time evolution of the specific an- 
gular momentum at the inner radius. The results 
are shown for the axisymmetric reference model R 
(triangles) and for the model A32 (squares, show- 
ing the range of values for different (^-directions). 
The filled symbols and solid lines denote the 9- 
averaged results, and the open symbols are for at 
the equatorial plane. 
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time [t orb (Rj] 

Fig. 12. — Time evolution of the mass accretion 
rate through the inner boundary (M; n ), in units of 
the Bondi rate (Mb). The non axisymmetric mod- 
els are: A (thick solid line), B (short dashed line), 
C (long dashed line), D (dotted line), E (dot- 
dashed line) and reference axisymmetric model R 
(thin solid line). 
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Fig. 13. — Time evolution of the angular momen- 
tum flux through the inner boundary, in units of 
the critical angular momentum times the Bondi 
accretion rate. The non axisymmetric models are: 
A (thick solid line), B (short dashed line), C (long 
dashed line) , D (dotted line) , E (dot-dashed line) 
and the reference axisymmetric model R (thin 
solid line). 



20 




Fig. 14.— The tilt angle (3 (Eq. [7) as a function 
of radius, for various times in the late phase of the 
evolution: t =0.18 (solid line), 0.21 (dotted line), 
0.25 (short dashed line), 0.29 (long dashed line), 
0.32 (dotted-short dashed line) and 0.34 (dotted- 
long dashed line). The initial tilt was zero (ini- 
tially axisymmetric model). 



Fig. 15. — The cumulative twist angle 7 (Eq. |8| as 
a function of time, for the late phase of the evolu- 
tion in the initially axisymmetric model. The solid 
line shows the twist averaged over innermost radii, 
from 1.5 x 10~ 3 to 5 x 10~ 2 Rb, while the dashed 
line shows the twist averaged over the whole disk, 
up to 1.2 i?B- 
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Fig. 16. — Density isosurfaces for five time snap- 
shots: 0.018, 0.09, 0.16, 0.23 and 0.29 t orh (R B ), for 
one of the 'broken torus' models, Bqq (see Tab. [TJ. 
The colors mark the surfaces of p = 500, 1250 and 
2000 poo . The boxes are scaled to the density iso- 
contour p — 500 p,^ and show the central region 
of the radius about 0.004 Rb- 
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